{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run a Gene Ontology Enrichment Analysis (GOEA)\n",
    "We use data from a 2014 Nature paper:    \n",
    "[Computational analysis of cell-to-cell heterogeneity\n",
    "in single-cell RNA-sequencing data reveals hidden \n",
    "subpopulations of cells\n",
    "](http://www.nature.com/nbt/journal/v33/n2/full/nbt.3102.html#methods)\n",
    "\n",
    "Note: you must have the Python package, **xlrd**, installed to run this example. \n",
    "\n",
    "Note: To create plots, you must have:\n",
    "  * Python packages: **pyparsing**, **pydot**\n",
    "  * [Graphviz](http://www.graphviz.org/) loaded and your PATH environmental variable pointing to the Graphviz bin directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Download Ontologies and Associations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1a. Download Ontologies, if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: go-basic.obo\n"
     ]
    }
   ],
   "source": [
    "# Get http://geneontology.org/ontology/go-basic.obo\n",
    "from goatools.base import download_go_basic_obo\n",
    "obo_fname = download_go_basic_obo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1b. Download Associations, if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: gene2go\n"
     ]
    }
   ],
   "source": [
    "# Get ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2go.gz\n",
    "from goatools.base import download_ncbi_associations\n",
    "gene2go = download_ncbi_associations()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load Ontologies, Associations and Background gene set "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2a. Load Ontologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "load obo file go-basic.obo\n",
      "46518"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "go-basic.obo: format-version(1.2) data-version(releases/2016-04-27)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " nodes imported\n"
     ]
    }
   ],
   "source": [
    "from goatools.obo_parser import GODag\n",
    "\n",
    "obodag = GODag(\"go-basic.obo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Load Associations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19,030 annotated mouse genes\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from goatools.associations import read_ncbi_gene2go\n",
    "\n",
    "geneid2gos_mouse = read_ncbi_gene2go(\"gene2go\", taxids=[10090])\n",
    "\n",
    "print(\"{N:,} annotated mouse genes\".format(N=len(geneid2gos_mouse)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c. Load Background gene set\n",
    "In this example, the background is all mouse protein-codinge genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from goatools.test_data.genes_NCBI_10090_ProteinCoding import GeneID2nt as GeneID2nt_mus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Initialize a GOEA object\n",
    "The GOEA object holds the Ontologies, Associations, and background.    \n",
    "Numerous studies can then be run withough needing to re-load the above items.    \n",
    "In this case, we only run one GOEA.    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fisher module not installed.  Falling back on scipy.stats.fisher_exact\n",
      "18,437 out of 28,212 population items found in association\n"
     ]
    }
   ],
   "source": [
    "from goatools.go_enrichment import GOEnrichmentStudy\n",
    "\n",
    "goeaobj = GOEnrichmentStudy(\n",
    "        GeneID2nt_mus.keys(), # List of mouse protein-coding genes\n",
    "        geneid2gos_mouse, # geneid/GO associations\n",
    "        obodag, # Ontologies\n",
    "        propagate_counts = False,\n",
    "        alpha = 0.05, # default significance cut-off\n",
    "        methods = ['fdr_bh']) # defult multipletest correction method\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Read study genes\n",
    "~400 genes from the Nature paper supplemental table 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Data will be stored in this variable\n",
    "import os\n",
    "geneid2symbol = {}\n",
    "# Get xlsx filename where data is stored\n",
    "ROOT = os.path.dirname(os.getcwd()) # go up 1 level from current working directory\n",
    "din_xlsx = os.path.join(ROOT, \"tests/data/nbt_3102/nbt.3102-S4_GeneIDs.xlsx\")\n",
    "# Read data\n",
    "if os.path.isfile(din_xlsx):  \n",
    "    import xlrd\n",
    "    book = xlrd.open_workbook(din_xlsx)\n",
    "    pg = book.sheet_by_index(0)\n",
    "    for r in range(pg.nrows):\n",
    "        symbol, geneid, pval = [pg.cell_value(r, c) for c in range(pg.ncols)]\n",
    "        if geneid:\n",
    "            geneid2symbol[int(geneid)] = symbol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Run Gene Ontology Enrichment Analysis (GOEA)\n",
    "You may choose to keep all results or just the significant results. In this example, we choose to keep only the significant results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating uncorrected p-values using fisher_scipy_stats\n",
      "   376 out of    400 study items found in association\n",
      "Running multitest correction: statsmodels fdr_bh\n",
      "16,953 GO terms are associated with 376 of 400 study items in a population of 28,212\n"
     ]
    }
   ],
   "source": [
    "# 'p_' means \"pvalue\". 'fdr_bh' is the multipletest method we are currently using.\n",
    "geneids_study = geneid2symbol.keys()\n",
    "goea_results_all = goeaobj.run_study(geneids_study)\n",
    "goea_results_sig = [r for r in goea_results_all if r.p_fdr_bh < 0.05]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Write results to an Excel file and to a text file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    164 items WROTE: nbt3102.xlsx\n",
      "    164 items WROTE: nbt3102.txt\n"
     ]
    }
   ],
   "source": [
    "goeaobj.wr_xlsx(\"nbt3102.xlsx\", goea_results_sig)\n",
    "goeaobj.wr_txt(\"nbt3102.txt\", goea_results_sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Plot all significant GO terms\n",
    "Plotting all significant GO terms produces a messy spaghetti plot. Such a plot can be useful sometimes because you can open it and zoom and scroll around. But sometimes it is just too messy to be of use.\n",
    "\n",
    "The **\"{NS}\"** in **\"nbt3102_{NS}.png\"** indicates that you will see three plots, one for \"biological_process\"(BP), \"molecular_function\"(MF), and \"cellular_component\"(CC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  WROTE: nbt3102_CC.png\n",
      "  WROTE: nbt3102_MF.png\n",
      "dot: graph is too large for cairo-renderer bitmaps. Scaling by 0.591505 to fit\n",
      "\n",
      "  WROTE: nbt3102_BP.png\n"
     ]
    }
   ],
   "source": [
    "from goatools.godag_plot import plot_gos, plot_results, plot_goid2goobj\n",
    "\n",
    "plot_results(\"nbt3102_{NS}.png\", goea_results_sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7a. These plots are likely to messy\n",
    "The *Cellular Component* plot is the smallest plot...\n",
    "![BIG CC PLOT](images/nbt3102_CC.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7b. So make a smaller sub-plot\n",
    "This plot contains GOEA results:\n",
    "  * GO terms colored by P-value:\n",
    "    * pval < 0.005 (light red)\n",
    "    * pval < 0.01 (light orange)\n",
    "    * pval < 0.05 (yellow)\n",
    "    * pval > 0.05 (grey) Study terms that are not statistically significant\n",
    "  * GO terms with study gene counts printed. e.g., \"32 genes\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  WROTE: nbt3102_MF_RNA_genecnt.png\n"
     ]
    }
   ],
   "source": [
    "# Plot subset starting from these significant GO terms\n",
    "goid_subset = [\n",
    "    'GO:0003723', # MF D04 RNA binding (32 genes)\n",
    "    'GO:0044822', # MF D05 poly(A) RNA binding (86 genes)\n",
    "    'GO:0003729', # MF D06 mRNA binding (11 genes)\n",
    "    'GO:0019843', # MF D05 rRNA binding (6 genes)\n",
    "    'GO:0003746', # MF D06 translation elongation factor activity (5 genes)\n",
    "]\n",
    "plot_gos(\"nbt3102_MF_RNA_genecnt.png\", \n",
    "    goid_subset, # Source GO ids\n",
    "    obodag, \n",
    "    goea_results=goea_results_all) # Use pvals for coloring\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RNA subplot](images/nbt3102_MF_RNA_genecnt.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7c. Add study gene Symbols to plot\n",
    "e.g., *11 genes: Calr, Eef1a1, Pabpc1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  WROTE: nbt3102_MF_RNA_Symbols.png\n"
     ]
    }
   ],
   "source": [
    "plot_gos(\"nbt3102_MF_RNA_Symbols.png\", \n",
    "    goid_subset, # Source GO ids\n",
    "    obodag,\n",
    "    goea_results=goea_results_all, # use pvals for coloring\n",
    "    # We can further configure the plot...\n",
    "    id2symbol=geneid2symbol, # Print study gene Symbols, not Entrez GeneIDs\n",
    "    study_items=6, # Only only 6 gene Symbols max on GO terms\n",
    "    items_p_line=3, # Print 3 genes per line\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RNA subplot](images/nbt3102_MF_RNA_Symbols.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2016, DV Klopfenstein, H Tang. All rights reserved."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
